Dynamics of a System of Two Coupled Oscillators Which Are Driven by a Third Oscillator

Analytical and numerical methods are applied to a pair of coupled nonidentical phase-only oscillators, where each is driven by the same independent third oscillator. The presence of numerous bifurcation curves defines parameter regions with 2, 4, or 6 solutions corresponding to phase locking. In all cases, only one solution is stable. Elsewhere, phase locking to the driver does not occur, but the average frequencies of the drifting oscillators are in the ratio of m:n. These behaviors are shown analytically to exist in the case of no coupling, and are identified using numerical integration when coupling is included. INTRODUCTION Recent experiments in optical laser MEMs have involved models of two coupled oscillators, each of which is being driven by a common harmonic forcer in the form of light [1]. Various steady states have been observed, including complete synchronization, in which both oscillators operate at the same frequency as the forcer, and partial synchronization, in which only one of the oscillators operates at the forcing frequency. Other possible steady states may exist, for example where the two oscillators are mutually synchronized ∗Address all correspondence to this author. but operate at a different frequency (or related frequencies) than that of the forcer (“relative locking”). Additionally, the oscillators may operate at frequencies unrelated to each other or to the forcing frequency (“drift”). The question of which of these various steady states is achieved will depend upon both the frequencies of the individual uncoupled oscillators relative to the forcing frequency, as well as upon the nature and strength of the forcing and of the coupling between the two oscillators. Related studies have been done for other variants of the three coupled oscillator problem. Mendelowitz et al. [2] discussed a case with one-way coupling between the oscillators in a loop; this system resulted in two steady states due to choice of direction around the loop. Baesens et al. [3] studied the general three–oscillator system (all coupling patterns considered), provided the coupling was not too strong, by means of maps of the two–torus. Cohen et al. [4] modeled segments of neural networks as coupled limit cycle oscillators and discussed the special case of two coupled phase–only oscillators as described by the following: θ̇1 = ω1 +α sin(θ2−θ1) (1) θ̇2 = ω2 +α sin(θ1−θ2) (2) 1 Copyright c © 2014 by ASME Defining a new variable ψ = θ2− θ1, being the phase difference between the two oscillators, allows the state of the system to be consolidated into a single equation: ψ̇ = ω2−ω1−2α sinψ (3) This is solved for an equilibrium point (constant ψ) which represents phase locking, i.e. the two oscillators traveling at the same frequency with some constant separation. This also gives us a constraint on the parameters which allow for phase locking to occur. If no equilibrium point exists, the oscillators will drift relative to each other; while they are affected by each other’s phase, the coupling is not strong enough to compensate for the frequency difference. sinψ∗ = ω2−ω1 2α (4) ∣∣∣∣ω2−ω1 2α ∣∣∣∣≤ 1 (5) Under the constraint, there are two possible equilibria ψ∗ and (π−ψ∗) within the domain, though one of them is unstable given that dψ̇ dψ =−2α cos(ψ) =−2α(−cos(π−ψ∗)) So if one of them is stable (dψ̇/dψ < 0), the other must be unstable (and vice versa). Plugging the equilibrium back into the original equations we find the frequency at which the oscillators end up traveling; this is a “compromise” between their respective frequencies. θ̇1 = ω1 +α ( ω2−ω1 2α ) = ω1 +ω2 2 (6) Since the coupling strength is the same in each direction, the resultant frequency is an average of the two frequencies with equal weight; with different coupling strengths this would become a weighted average. Keith and Rand [5] added coupling terms of the form α2 sin(θ1−2θ2) to this model and correspondingly found 2:1 locking as well as 1:1 locking. MODEL We design our model, as an extension of the two– oscillator model, to include a pair of coupled phase-only oscillators with a third forcing oscillator, as follows: θ̇1 = ω1 +α sin(θ2−θ1)−β sin(θ1−θ3) (7) θ̇2 = ω2 +α sin(θ1−θ2)−β sin(θ2−θ3) (8) θ̇3 = ω3 (9) This system can be related back to previous work by careful selection of parameters. Note that the β = 0 case reduces the system to two coupled oscillators without forcing, while α = 0 gives a pair of uncoupled forced oscillators. It is now useful to shift to a coordinate system based off of the angle of the forcing oscillator, since its frequency is constant. Let φ1 = θ1− θ3 and φ2 = θ2− θ3, with similar relations Ω1 = ω1−ω3 and Ω2 = ω2−ω3 between the frequencies. The forcing oscillator’s equation of motion can thus be dropped. φ̇1 = Ω1 +α sin(φ2−φ1)−β sinφ1 (10) φ̇2 = Ω2 +α sin(φ1−φ2)−β sinφ2 (11) A nondimensionalization procedure, scaling time with respect to Ω1, allows for Ω1 = 1 to be assumed without loss of generality. We note that the φi now represent phase differences between the paired oscillators and the driver. Thus, if a φ̇i = 0, the corresponding θi is defined to be locked to the driver. Equilibrium points of equations (10) and (11) then represent full locking of the system. Partial and total drift are more difficult to recognize and handle analytically, and will be discussed later. FULL LOCKING We begin by solving the differential equations for equilibria directly, so as to find the regions of parameter space for which the system locks to the driver. The equilibria satisfy the equations: 0 = 1+α sin(φ2−φ1)−β sinφ1 (12) 0 = Ω2 +α sin(φ1−φ2)−β sinφ2 (13) Trigonometrically expanding equation (12) and solving for cos φ1: cos φ1 = α sin φ1 cos φ2 +β sin φ1−1 α sin φ2 (14) We square this equation, rearrange it, and use sin2 θ + cos2 θ = 1 to replace most cosine terms: α2 (1− sin2 φ1) sin2 φ2 −α2 sin2 φ1(1− sin2 φ2) +(2α sin φ1−2αβ sin2 φ1) cos φ2 −β2 sin2 φ1 +2β sin φ1−1 = 0 (15) Repeating the process by solving for cos φ2, we obtain an equation in terms of only sines: −[4α2β2 sin4 φ1−8αβ sin3 φ1 2 Copyright c © 2014 by ASME +(4α2−2α2β2−2α4) sin2 φ1 +4α2β sin φ1−2α] sin2 φ2 −α4 sin4 φ2− (β4−2α2β2 +α4) sin4 φ1 −(4α2β−4β3) sin3 φ1− (6β2−2α2) sin2 φ1 +4β sin φ1−1 = 0 (16) Returning to equations (12) and (13), we add them and solve for sinφ2 in terms of sinφ1: sin φ2 = 1+Ω2 β − sin φ1 (17) This is now plugged into equation (16) to eliminate φ2 and obtain an polynomial of degree six in s = sin φ1, dependent on the various parameters. −4α2β6s6 +8αβΩ2s +16αβs−4αβΩ2s −24αβΩ2s−βs +4α2β6s4−24α2β4s4 +8αβΩ2s −4αβΩ2s +24αβΩ2s +4β7s3 −12α2β5s3 +16α2β3s3 +2αβΩ2s−4αβΩ2s −4αβΩ2s +12αβΩ2s−8αβΩ2s −8αβΩ2s−6βs +14α2β4s2−4α4β2s2 −4α2β2s2 +4αβΩ2s−4αβΩ2s+12αβΩ2s −12αβΩ2s+12αβΩ2s+4βs−8αβs +4αβs−αΩ2−4αΩ2 +2αβΩ2−6αΩ2 +4αβΩ2−4αΩ2−β +2α2β2−α4 = 0 (18) The roots of this polynomial give values of s = sin φ1 for a given set of parameter values; degree six implies that there will be up to six roots in s, although a single root in s may correspond to more than one root in φ1 due to the multivalued nature of sine. Each φ1 has a corresponding φ2 value as defined by equation (17). To avoid extraneous roots, each φi pair should be confirmed in equations (12) and (13). In order to distinguish changes in the number of real equilibria, we look for double roots of this polynomial such that two (or more) of the equilibria are coalescing in a single location. Setting the polynomial and its first derivative in s equal to zero and using Maxima to eliminate s results in a single equation with 142 terms in α, β, and Ω2 which describes the location of bifurcations. 48Ω2α β−360Ω2αβ−32Ω2αβ−87Ω2α +64Ω10 2 α 6 +320Ω2α −128Ω2αβ −1308Ω2αβ−112Ω2αβ−160Ω2α +616Ω2α 6 +16Ω2α 4 +12Ω2α 2β8 +22Ω2α 4β6 −16Ω2αβ +410Ω2αβ +528Ω2α +64Ω2α +6Ω2α 4β4 +8Ω2α 2β4 +140Ω2α 8β2 −1720Ω2αβ−224Ω2αβ +256Ω2α −52Ω2αβ +324Ω2αβ +100Ω2αβ +796Ω2α 8 +88Ω2α 6 +96Ω2α 4 +1952Ω2α 6β4 +448Ω2α β−48Ω2αβ−1240Ω2αβ −996Ω2αβ−288Ω2αβ +Ω2β −34Ω2αβ +1536Ω2α +3232Ω2α −160Ω2α +64Ω2α−2Ω2β−189Ω2αβ +54Ω2α 2β8 +Ω2β −480Ω2αβ +94Ω2α β−40Ω2αβ +960Ω2αβ +1526Ω2α 6β4 +404Ω2α 4β4 +8Ω2α 2β4 −1024Ω2αβ−6284Ω2αβ−448Ω2αβ −224Ω2αβ +3840Ω2α +4726Ω2α +88Ω2α 6 +108Ω2α β−152Ω2αβ −208Ω2αβ +16Ω2α−752Ω2αβ +100Ω2α 2β6 +2560Ω2α 8β4 +48Ω2α 6β6 −96Ω2αβ +448Ω2αβ−4096Ω2αβ −9808Ω2αβ−996Ω2αβ−112Ω2αβ +5120Ω2α 10 +3232Ω2α 8 +528Ω2α 6 −2Ω2β +30Ω2αβ +68Ω2αβ −56Ω2αβ−2Ω2β−320Ω2αβ −166Ω2αβ +54Ω2αβ +512Ω2αβ +864Ω2α 6β6 +94Ω2α β−16Ω2αβ +3200Ω2α 8β4 +1526Ω2α 6β4 +6Ω2α 4β4 −6144Ω2αβ−6284Ω2αβ−1720Ω2αβ −56Ω2αβ +152Ω2αβ−32Ω2αβ +3840Ω2α 10 +4Ω2β 12 +796Ω2α 8 +616Ω2α 6 +108Ω2αβ−152Ω2αβ −52Ω2αβ +1024Ω2αβ +48Ω2αβ +324Ω2αβ +2560Ω2αβ +1952Ω2αβ −128Ω2αβ−832Ω2αβ−4096Ω2αβ −1240Ω2αβ +β−12αβ−1308Ω2αβ +1536Ω2α−160Ω2α +320Ω2α−2β +48α4β12 +30α2β12 +β12−64α6β10 +68α4β10 −34α2β10−320α6β8−189α4β8 +12α2β8 +512α8β6−480α6β6 +22α4β6 +960α8β4 +410α6β4 +48α4β4−1024α10β2 +140α8β2 −360α6β2 +256α10−87α8 +64α6 = 0 (19) This equation by itself is cumbersome to work with. We begin to interpret its results by choosing different values of Ω2 and plotting the resulting curves in the βα-plane (see Fig. 1). Each curve is the location of a double root of the original system, and represents a pair of equilibrium points being created or destroyed in a fold bifurcation. The combination of bifurcation curves leads to regions of 0–6 equilibria. Numerical analysis of the original differential equations with AUTO continuation software [6] both confirms the quantities of equilibria and calculates the eigenvalues of each point. Through these results, we find that only one equilibrium point (and therefore locking behavior) is stable; it occurs for any region where equilibria exist, i.e. for large 3 Copyright c © 2014 by ASME 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14